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Abstract 



We comment on the reweighting method for the study of finite density lattice 
QCD. We discuss the applicable parameter range of the reweighting method 
for models which have more than one simulation parameter. The applica- 
bility range is determined by the fluctuations of the modification factor of 
the Boltzmann weight. In some models having a first order phase transition, 
the fluctuations are minimized along the phase transition line if we assume 
that the pressure in the hot and the cold phase is balanced at the first order 
phase transition point. This suggests that the reweighting method with two 
parameters is applicable in a wide range for the purpose of tracing out the 
phase transition line in the parameter space. To confirm the usefulness of 
the reweighting method for 2 flavor QCD, the fluctuations of the reweight- 
ing factor are measured by numerical simulations for the cases of reweighting 
in the quark mass and chemical potential directions. The relation with the 
phase transition line is discussed. Moreover, the sign problem caused by the 
complex phase fluctuations is studied. 

11.15.Ha, 12.38.Gc, 12.38.Mh 
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I. INTRODUCTION 



The study of the phase structure of QCD at non-zero temperature T and non-zero quark 
chemical potential \x q is currently one of the most attractive topics in particle physics [1,2]. 
The heavy-ion collision experiments aiming to produce the quark-gluon plasma (QGP) are 
running at BNL and CERN, for which the interesting regime is rather low density. Moreover 
a new color superconductor phase is expected in the region of low temperature and high 
density. In the last few years, remarkable progress has been achieved in the numerical study 
by Monte-Carlo simulations of lattice QCD at low density. It was shown that the phase 
transition line, separating hadron phase and QGP phase, can be traced out from /i q — 
to finite fi q , and it was also possible to investigate the equation of state quantitatively at 
low density. The main difficulty of the study at non-zero baryon density is that the Monte- 
Carlo method is not applicable directly at finite density, since the fermion determinant is 
complex for non-zero ji q and configurations cannot be generated with the probability of the 
Boltzmann weight. The most popular technique to study at non-zero /i q is the reweighting 
method; performing simulations at Re(/x 9 ) = 0, and then modify the Boltzmann weight at the 
step of measurement of observables [3-6] . The Glasgow method [7] is one of the reweighting 
methods. A composite (Glasgow) reweighting method has recently been proposed by [8]. 
Another approach is analytic continuation from simulations at imaginary chemical potential 
[9-1 1] . Moreover, calculating coefficients of a Taylor expansion in terms of \i q is also a hopeful 
approach for the study at non-zero baryon density [4,12-14]. The studies by Taylor expansion 
or imaginary chemical potential require analyticity of physical quantities as functions of T 
and n q , while the reweighting method has a famous "sign problem". The sign problem 
is caused by complex phase fluctuations of the fermion determinant, which are measured 
explicitly in Ref. [15], and also Ref. [16] is a trial to avoid the sign problem. 

In this paper, we make comments on the reweighting method with respect to more than 
one simulation parameter, particularly, including a chemical potential. Because the fluc- 
tuation of the modification factor (reweighting factor) enlarges the statistical error, the 
applicable range of this method is determined by the fluctuation of the reweighting fac- 
tor. When the error becomes considerable in comparison with the expectation value, the 
reweighting method breaks down. An interesting possibility is given in the following case: 
when we change two parameters at the same time, the reweighting factors for two param- 
eters might cancel each other, then the error does not increase and also the expectation 
values of physical quantities do not change by this parameter shift. Therefore finding such a 
direction provides useful information for mapping out the value of physical quantities in the 
parameter space. As we will see below, there is such a direction in the parameter space and 
the knowledge of this property of the reweighting method makes the method more useful. 
Fodor and Katz [3] investigated the phase transition line for rather large fi q . This argument 
may explain why they could calculate f3 c for such large \i q . 

In the next section, we explain the reweighting method with multi-parameters. Then, 
in Sec. Ill, the case of SU(3) pure gauge theory on an anisotropic lattice is considered as 
the simplest example, and we proceed to full QCD with a first order phase transition in 
Sec. IV. For 2-flavor QCD, the reweighting method with respect to quark mass is discussed 
in Sec. V. The application to non-zero baryon density is discussed in Sec. VI. The problem 
of the complex measure is also considered in Sec. VII. Conclusions are given in Sec. VIII. 
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II. REWEIGHTING METHOD AND THE APPLICABILITY RANGE 



The reweighting method is based on the following identity: 

(0) {M = I DUO(detM(m,i,)) N <e- s ^ 

(Oe F e G ) (Oe AF e AG ) 

2 / (/3o,mo,fM)) _ \ / (/3o,mo,/jo) /, \ 

~~ (e F e G ) ~ (e AF e AG ) ' 

\ '(j3o,mo,no) \ ' (/3o, mo, no) 

Here M is the quark matrix, S g is the gauge action, Nf is the number of flavors (iVf/4 
for staggered type fermions instead of Nf), F = Nf [In det M(m, //) — In det M(m , /io)], 
G — (P — p )P,P = -dS g /df3, (3 = 6/g 2 , and AF[G] = - m and /i = // 9 a 

are the bare quark mass and chemical potential in a lattice unit, respectively. The expec- 
tation value (0)(p >m>fl -) can, in principle, be computed by simulation at (p , m , /i ) using 
this identity [17]. However, a problem of the reweighting method is that the fluctuation 
of e AF e AG = e F e G /(e^e^) enlarges the statistical error of the numerator and denom- 
inator of Eq.(l). The worst case, which is called "Sign problem", is that the sign of the 
reweighting factor changes frequently during Monte-Carlo steps, then the expectation values 
in Eq.(l) become vanishingly small in comparison with the error, and this method does not 
work. However, the fluctuations cause the f3, m and \i dependence of {0)^ >m ^). Otherwise, 
if e F e G does not fluctuate, e AF e AG = 1 and (0)(p >m , ifl ) does not change with parameter 
change. Roughly speaking, the difference of (O)^,™^) from (C , )( / 3 , m . 0)W ) increases as the 
magnitude of fluctuations of F and G increases and, if F and G have a correlation, the 
increase of the total fluctuation as a function of f3, m and ji is non-trivial. Therefore, it is 
important to discuss the correlation between e F and e G and to estimate the total fluctuation 
of the reweighting factor in the parameter space in order to estimate the applicability range 
of the reweighting method in the parameter space and how the system changes by parameter 
shifts. This is helpful information for the study of QCD thermodynamics. 



III. SU(3) GAUGE THEORY ON AN ANISOTROPIC LATTICE 

Let us start with the case of SU(3) pure gauge theory on an anisotropic lattice, having 
two different lattice spacings for space and time directions: a a and a T . As we will show for 
this case, there is a clear relation between the phase transition line and the direction which 
minimizes the fluctuation of the reweighting factor. The action is 

Sg = -Pa E P M ~ Pr £ P T {x\ (2) 

X X 

where P a (r) is spatial (temporal) plaquette. The SU(3) pure gauge theory has a first order 
phase transition. At the transition point (T c ), there exist two phases simultaneously. For the 
two phases to coexist, the pressure in these phases must be equal: Ap = p( hot ) — p( cold ) = Q, 
If we require Ap = 0, we find that the phase transition line in the parameter space of (P a , P T ) 
has to run in such a direction that the fluctuation of the reweighting factor is minimized 
when we perform a simulation on the phase transition point and apply the reweighting in 
the (P a ,P T ) plane. 
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A significant feature of a Monte-Carlo simulation at a first order phase transition point 
is the occurrence of flip-flops between configurations of hot and cold phases. If one write a 
histogram of the action density, i.e. the plaquette: P a and P T , there exist two peaks. The 
value of the action density changes sometimes from one near one peak to one near the other 
peak during Monte-Carlo steps [18,19]. The flip-flop is the most important fluctuation; 
in fact, the flip-flop makes the strong peak of susceptibilities of observables such as the 
plaquette or the Polyakov loop at the transition point. Also, the flip-flop implies strong 
correlations between P CT and P T because the values of P a and P T change simultaneously 
between the typical values of the two phases in the (P CT , P T ) plane. 

Here, we discuss the fluctuation of the reweighting factor when one performs a simulation 
at the transition point, (/3 CT o, Pro)- An expectation value of O at (/3 a , f3 T ) on an anisotropic 
lattice is calculated by 



(O) 



09a A) 



Oe 



-AS, 



(0*0, Pro)/ 



e ~AS a 



(3) 



where AS g = —Af3 a J2 X Pa{x) — A(5 T J2 X Pt(x), and A/3 CT ( T ) = /3 CT ( T ) — /3 CT (r)o- For simplification, 
we ignore the local fluctuation around the two peaks of the histogram of the action density 
and consider only the flip-flop between hot and cold phases, since it is the most important 
fluctuation at the first order phase transition point. Then, the fluctuation is estimated by 
the difference of the reweighting factor between hot and cold phases, up to first order, 



,(hot) 



(cold) 



3N site \Af3 a (P^ - Pi cold) ) + A/3 r (Pf ot ) - P> old ))| + • • • , (4) 



where phot(coid) and phot(coid) 

are the average values of the spatial and temporal plaquettes 
for configurations in the hot (cold) phase and iV S i te = N% x N T is the number of sites for an 
N% x N T lattice. Hence, along the line which has a slope 



Jj(hot) p(co\d) 



PoV 



(hot) 



PS 



(cold) ' 



(5) 



the fluctuation of the reweighting factor is canceled to leading order. 

On the other hand, since V = (Af ff o ff ) 3 and T = (A^" T a r ) _1 , pressure is defined by 



p = T 



dlnZ 



dV 



1 



Ni 



3N*N T ala T da a 
d(5 a 



d(3 Ti 



, . )a -({Pa) - (P,) ) + Oa-^m) - (Prh) 



(6) 
(7) 



where Pa(r) = (SA^ite) -1 J2 X Pa(r)(x), and (P a ( T ))o is the expectation value of the plaquette 
on a T = lattice for the normarization. 

By separating the configurations into those in hot and cold phases [18,19], the gap of 
pressure between hot and cold phases at T c is computed by 



p p( hot ) 
= iV r 4 (— 



V 



(cold) 



rp4 



"a 



p(cold)> 



p(cold) n 



(8) 



4 



Since the gap of pressure should vanish, Ap = 0, 



df^/dPr P r (hot) - P T (cold) ( ] 

Oa J da a pjhot) _ pMd) ■ W 

Moreover, because T c = (Af r a r ) _1 on the phase transition line, a T keeps constant with 
(NtTc)^ 1 along the transition line, i.e. 

Aa T = ^-A(3 a + ^A[3 T = 0, (10) 
o(3 c dj3 T 

when one changes parameters ((3 a ,(3 T ) — > (f3 a + Af3 a) f3 T + Af3 T ) along the phase transition 
line. Then the slope of the phase transition line (r t ) [20] is obtained by 



d/3 CT 

n = 



d(3. 

where we used an identify: 



da T I da T df3 a / d(3 T 



df3 T df3 a da a da^ 



( df3 a d/3 T \ i / da T da T \ 

d£o_ dPr I - da T day da T da a I _§a 3 _ daa_ I ■ \ LZ ) 
da T da T I d/3 T d/3a d(3 a df3 T V dPr d(3 a / 

Hence, the condition for Ap = becomes 

p(hot) p(cold) 

~^t) K^Mj = _r *- ( 13 ) 

Co Co 

This equation for d(5 c /d(5 T may correspond to the Clausius-Clapeyron equation in the 
(p,T) plane: dp/dT = AS/AV (S: entropy). In fact, Eq.(13) and the vanishing pressure 
gap, Ap = 0, are confirmed by calculating the slope of the transition line from the peak 
position of the Polyakov loop susceptibility obtained by numerical simulations in Ref. [20]. 
Historically, the non-zero gap of pressure at the transition point had been a problem for 
long time. The reason of Ap ^ was that the precise non-perturbative measurement 
of the anisotropy coefficients, a a (dfl a /da a ),a a (dfl T /da a ), etc., had been difficult, and the 
perturbative value [21] cannot be used at the phase transition point for iV T = 4 or 6. After 
the precise measurement of the anisotropy coefficients became possible, the problem of a 
non-zero pressure gap was solved. Also, determinations of these coefficients by another 
non-perturbative method have been done in Refs. [22-24]. 

From Eq.(5) and Eq.(13), we find that the direction which minimizes the fluctuation 
of the reweighting factor must be the same as the direction of the phase transition line in 
the (f3 a ,(3 T ) plane, if pressure in hot and cold phases are balanced at T c , Ap = 0. Here, 
in practice, we estimate the fluctuation of e~ ASs > by a numerical simulation. We compute 
the standard deviation of the reweighting factor using the data obtained in Ref. [19]. The 
lattice size is 24 2 x 36 x 4. The data is generated by the standard Wilson gauge action 
with (3 a = (3 a = 5.6925, that is just on the transition line, f3 c = 5.69245(23) for A^ T = 4 
at a a = a T . Ellipses in Fig.l are the contour lines of the standard deviation normalized 
by the mean value: y / ((e _AS 's) 2 ) — (e~ AS 3) 2 / (e~ AS9 ), and we write this value in Fig. 1. We 
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also denote the phase transition line, obtained by the measurement of the Polyakov loop 
susceptibility assuming that the peak position of the susceptibility is the phase transition 
point [20], by a bold line, and dashed lines are the upper bound and lower bound. We 
find that the phase transition line and the line which minimizes fluctuations are consistent. 
This result also means that the reweighting method is applicable in a wide range for the 
determination of the phase transition line in the parameter space of SU(3) gauge theory on 
an anisotropic lattice, since the increase of the statistical error caused by the reweighting is 
small along the transition line. 

We note that, from Eq.(13), AS^ hot) and A^ cold) in Eq.(4) are equal under the change 
along the phase transition line, hence the fluctuations are canceled in every order of A/3 CT ( r ). 
For SU(3) pure gauge theory on an anisotropic lattice, the system is independent of a a /a T 
in a physical unit, hence the system does not change along the transition line except for the 
volume, V = (N^a^) 3 , if N a is finite. 1 Because physical quantities does not change without 
the fluctuation of reweighting factor, this result is quite natural. 



IV. FULL QCD WITH A FIRST ORDER PHASE TRANSITION 



Next, we extend this discussion for the case of full QCD with a first order phase transition 
such as 3 flavor QCD near the chiral limit. The reweighting method is applied in the 
parameter space of (m,(3). We consider Helmholtz free energy density / = — T\nZ/V for 
a canonical ensemble which is equal to minus pressure, p = —f, for a large homogeneous 
system. Under a parameter change from (m, (5) to (m + Am, (5 + A/5), the variation of the 
free energy is given, up to the first order, by 



f_ 



(m+Am,/3+A/3) 



rp4 



-Nt 



(m,/3) 



«0i> - (Oi>o)Am + 6((P) - (P) )A(3 



+ 



(14) 



where Q\ = N~^ c N { <9(ln det M) /dm and P = — (6A? sitc ) _1 (9S' £ ,/<9/3. For the normalization at 
T = 0, we subtract the zero temperature contribution, (P) and (Qi)o- Here, we should 
notice that the first derivatives of the free energy are discontinuous at the phase transition 
line, hence we cannot estimate the difference of the free energy beyond the transition line 
by this equation. 

We assume that the gap of the pressure is zero in the entire parameter space (m, (5). We 
change m and (3 along the phase transition line starting at two points: just above and just 
below (3 C , without crossing the transition line. Then the change of the free energy must be 
the same for both these cases, since a pressure gap is not generated under this variation, i.e. 
A(Ap) = 0. Hence, up to first order of Am and A/5, 



(hot) 



(m+Am,/3+A/3) 



f_ 



(hot) 



(m,/3), 



(cold) 



(m+Am,/3+A/3) 



-Nt 



(gS hot) - g! cold) )A™ + 6(p (hot) - P^)Af3 



f_ 

+ ■ 



(cold) N 



(15) 



1 In fact, as expected from the finite size scaling, the peak height of the Polyakov loop susceptibility 
increases as a a increases [20]. 
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is required on the first order phase transition line. From this equation, we obtain a similar 
relation as Eq.(13): 



g(hot) _ QMd) 



d(3 



(16) 



6(P(hot) _ p(cold)) 




On the other hand, the change of the reweighting factor under flip-flop is 



AF( hot ) AG< hot ) AF( cold ) AG( cold ) 



iVsite [(QS hot) - g! cold) )A^ + 6(P (hot) - P (cold) )A/?] + • • • . 



(17) 



If we ignore the local fluctuation around the peaks of the distribution of P and Q\, again, 
the direction for which the fluctuation is canceled is 



This is the same direction as the phase transition line. Therefore, the fluctuation of the 
reweighting factor along the phase transition line remains small, i.e. the statistical error 
does not increase so much. 

We obtained the same result as for the pure gauge theory on an anisotropic lattice, 
and this argument seems to be quite general for models with a first order phase transition, 
including models with chemical potential. However, there is a difference. Under the change 
of a a /a T , any physics does not change along the T c line, but physical quantities, in general, 
depend on the quark mass. Although the dependence on m might be much smaller than the 
dependence on T/T c , if the fluctuation of the reweighting factor is completely canceled, any 
m-dependence is not obtained. In this discussion, we ignored the local fluctuation around 
the peaks in the hot and cold phase respectively, but the local fluctuation may play an 
important role for the m-dependence. Also, the sign problem for nonzero baryon density is 
caused by complex phase fluctuations of the reweighting factor (see Sec. VII.), that is by the 
local fluctuations. Hence the local fluctuation may, in particular, be important at non-zero 
baryon density. 



As we have seen in the previous two sections, the multi-parameter reweighting seems to 
be efficient to trace out the phase transition line in a wide range of the parameter space. One 
of the most interesting applications is finding the (pseudo-) critical line (f3 c ) in the (m, (5) 
plane for 2 flavor or 2+1 flavor QCD. The phase transition for 2 flavor QCD at finite quark 
mass is expected to be crossover, which is not related to any singularity in thermodynamic 
observables, and that for 3 flavor QCD is crossover for quark masses larger than a critical 
quark mass, and is of first order for light quarks. The precise measurement of the (pseudo-) 
critical line is required for the extrapolation to the physical quark masses and for the study of 
universality class, e.g. to investigate for the 2 flavor case whether the chiral phase transition 
at finite temperature is in the same universality class as the 3-dimensional 0(4) spin model 
or not. 



d/3 Qr } - Q [ l 



(18) 



dm 6(P( hot ) - P( cold )) 



V. QUARK MASS REWEIGHTING FOR 2 FLAVOR QCD 
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In Ref. [4], we applied the reweighting method in the (m, /3) plane for 2 flavor QCD, 
and calculated the slope of the transition line, d(3 c /dm, where the reweighting factor with 
respect to quark mass was expanded into a power series and higher order terms which does 
not affect the calculation of the slope were neglected. The results for d/3 c /dm compared well 
with the data of (3 c (m) obtained by direct calculations, without applying the reweighting 
method, demonstrating the reliability of a reweighting in a parameter of the fermion action. 
In this section, we discuss the relation between the fluctuation of the reweighting factor and 
the phase transition line in the (m, (3) plane for 2 flavor QCD at finite quark mass, i.e. at 
the crossover transition, by measuring the fluctuation in numerical simulations. 

For the presence of a direction for which the two reweighting factors from the gauge and 
the fermion action cancel, a correlation between these reweighting factors during Monte- 
Carlo steps is required. We estimate the correlation between these reweighting factors using 
the configurations in Ref. [4]. A combination of the Symanzik improved gauge action and 
2 flavors of the p-4 improved staggered fermion action is employed [25]. The parameters 
are m = 0.1, p = {3.64, 3.645, 3.65, 3.655, 3.66, 3.665 and 3.67}. The lattice size is 16 3 x 4. 
7800-58000 trajectories are used for measurements at each (3. The details are given in 
Ref. [4]. 2 

In the vicinity of the simulation point, the correlation of the reweighting factors can be 
approximated by 

(e F e G ) - (e F )(e G ) = (Ae F Ae G ) » (A<2iAP)(m - m )(P - ft) + • • • , (19) 

where P = -dS g /d(3, F = £^° =1 Q n (m - m ) n , and we denote AX = X - (X) for X = 
{P, Q n , ■ ■ •}. The Q n are obtained by 

Q 1 = (Ar f /4)trM-\ Q 2 = -{Nt/SMM^M- 1 ), • • • , (20) 

for standard staggered fermions and also for p4-improved staggered fermions. We calculate 
the value of (AQiAP) = (QiP) — (Q\)(P). The random noise method is used for the 
calculation of Q n . The results for (AQiAP) are listed in Table I. We find strong correlation 
between the gauge and fermion parts of the reweighting factor. 

Then we compute the total fluctuation of the reweighting factor as a function of m and 
P near the simulation point. Up to second order in (3 — p and m — m , the square of the 
standard deviation is written as 

([A(e F e G )] 2 ) » {(AQ 1 ) 2 )(m - m ) 2 + 2(AQ 1 AP)(m - m )((3 - (3 ) 

+ ((AP) 2 )( / 3-A ) ) 2 + ---. (21) 

If we approximate in this form, lines of constant fluctuation (standard deviation) in the 
(m, P) plane form ellipses. We also compute ((AQi) 2 ) and ((AP) 2 ), which are written in 
Table I. The values at P c = 3.6492(22) are interpolated by applying the reweighting method 
for P direction combining the data at seven simulation points [17]. The lines of constant 



2 The coefficient c F of the knight's move hopping term was incorrectly reported to be 1/96 in 
Ref. [4]; its correct value is 1/48. 
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fluctuation are drown in Fig. 2. Numbers in this figure are the squares of the standard 
deviation divided by N s[tc . It is found that these ellipses spread over one direction and the 
increase of the fluctuation is small along the direction. We also show the slope of the phase 
transition line by two lines: upper bound and lower bound of the derivative of f3 c with respect 
to m obtained by measuring the peak position of the chiral susceptibility, d(3 c /dm = 1.05(14) 
for itlq = 0.1 [4]. We see that the directions of small fluctuations and of the phase transition 
are roughly the same. Since the fluctuation enlarges the statistical error of an observable, 
this figure can be also regarded as a map indicating the increase of the statistical error due 
to the reweighting. Therefore, we understand that the reweighting method can be applied 
in a wide range of parameters along the phase transition line if one performs simulations at 
the phase transition point. 

Moreover, it might be also important that these two directions are not exactly the same 
because, if the fluctuation is completely canceled along the transition line, no quantity can 
change, but the system should change as a function of quark mass even on the transition 
line, e.g. the chiral susceptibility should become larger as m decreases. 

VI. CHEMICAL POTENTIAL REWEIGHTING FOR 2-FLAVOR QCD 
A. Correlation among the reweighting factors 

Next, let us discuss the reweighting method for no n- vanishing chemical potential. The 
reweighting method is really important for the study of finite density QCD since direct 
simulations are not possible for non-zero baryon density at present. However, the complex 
measure problem (sign problem) is known to be a difficult problem. The reweighting factor 
for non-zero \x is complex. If the complex phase fluctuates rapidly and the reweighting 
factor changes sign frequently, the expectation values in Eq.(l) become smaller than the 
error. Then the reweighting method breaks down. Therefore it is important to investigate 
the reweighting factor, including the complex phase, in practical simulations. 

First of all, we separate the fermion reweighting factor e F into an amplitude \e F \ and 
a phase factor e , and investigate the correlation among \e F \, e 10 and the gauge part e G , 
where e G is real. As is shown in Ref. [4], the phase factor and the amplitude can be written 
as the odd and even terms of the Taylor expansion of In det M, respectively, since the odd 
terms are purely imaginary and the even terms are real at /i — 0. Denoting F = X^li RnfJ> n , 

oo oo 

\e F \ = exp{]T Re(R 2n )fi 2n } and e ie = exp{* ]T lm(R 2n ^) fi 2 ^ 1 } . (22) 

n=l n=l 

We study the correlation among these factors in the vicinity of the simulation point (/3q, Ho = 
0). Up to 0{(5 — (3o,fi 2 ), the reweighting factor is 

e w \e F \e G &1 + R lf ji+ (P 2 /2)/i 2 + R 2f i 2 + P{(3 - #>) + • • • . (23) 

We compute the correlations, (A(i? 2 /2)AP), {AR 2 AP), and (A(P 2 /2)AP 2 ) at (i = 0, which 
correspond to the correlations of (e l6 ,e G ), (|e F |,e G ) and (e' t9 ,\e F \), respectively. AX = 
X - (X) for X = {P, R n , ■ ■ •}. Here, (Ai^AP) and (ARiAR 2 ) are zero at n = because 
i?i is purely imaginary. The R n are obtained by 
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Ri 
R2 
R3 



— tr M — — 
4 \ 0u 



7V f 9 In detM 

iVf 1 d 2 In det M _ N { 1 

iV f 1 9 3 (In detM) _ JV f 1 
~4~3! 9^ ~ T3! 



tr M 



_ x d 2 M\ _ ( x dM x dM" 
dfi 2 J \ " dfj, dfi 



(24) 
(25) 



tr ( M "3 tr (M -'^m 



+2tr M — — M — — M 



(26) 



for staggered type fermions. Details of the calculation are given in Ref. [4]. 

We use the configurations in Ref. [4], again, generated by the Nf — 2 j94-improved 
action on a 16 3 x 4 lattice. We generated 20000-40000 trajectories for m = 0.1, (3 — 
{3.64,3.65,3.66 and 3.67}. The results are summarized in Table II. We find that the 
correlation between \e F \ and e G is very strong in comparison with the other correlations, 
which means that the contribution to an observable can be separated into two independent 
parts: from e 10 , and from a combination of \e F \ x e G . 

To make the meaning of this result clearer, we consider the following partition function, 
introducing two different //; ji Q and fj, e , 



Z = J J Df/e i?1 ^ +i?s ^ + -e i?2 ^ +R4 ^ + -(detM|^ =0 ) 7Vi/4 e" 



Then, at ji — 0, 



(ARjAP) = ((Ai?i) 2 AP) = 
2(AR 2 AP) = 



d 3 \nZ 
d 3 lnZ 



site 



site" 



d(x q a 2 ~ 4xiv« 2 
dp 

dp ' 



(27) 

(28) 
(29) 



where Xq an d Xiv are the quark number susceptibility and iso-vector quark number suscep- 
tibility [26]: 



J^2 



d 



+ 



d 



n u + n d 



d^jT) d(n A /T)) T 3 

4xiv = / g d \ n u -n d 

T 2 ~ \d(n u /T) dQii/T)) T 3 ■ 

We choose the same chemical potential for up and down quarks: ji u = ji A = /i q . 



n 



u(d) 



(30) 

(31) 
is the 



number density for up (down) quark: n u ^/T 3 = <9(p/T 4 )/<9(/i u ( d )/T). Quark and baryon 
number susceptibilities are related by xb = dris/d^B = 3~ 2 x q - If we impose a chemical 
potential with opposite sign for up and down quarks: fi u = —fid = l^iv/^, called "iso-vector 
chemical potential" , the Monte-Carlo method is applicable since the measure is not complex 
[27,28]. For this model, the iso-vector quark number susceptibility xiv in Eq.(31) is the 
quark number susceptibility instead of Eq.(30). 
The result in Table II means that 



d 3 In Z d 3 In Z 

dfi 2 Q dp < dfi 2 e dp ' 



(32) 
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i.e. /x in the phase factor (/x G ) does not contribute to the /^-dependence of Z near /x = 0, 
hence /x in the amplitude (/x e ) is more important for the determination of ft by measuring 
the ftdependence of thermodynamic quantities. Moreover, these correlations have a relation 
with the slope of x q and 4xrv m terms of ft Since Xq ~~ 4%rv is known to be small at /x — 
[26] , 3 this result may not change even for small quark mass. Also, the result in Ref. [30] 
suggests that the effect of the phase factor, i.e. of /x G , on physical quantities is small. 

Iso-vector chemical potential Furthermore, we discuss the model with iso-vector chemi- 
cal potential. In Ref. [4], we discussed the difference in the curvature of the phase transition 
to that of the usual chemical potential. Because we expect that at T = pion condensation 
happens around /x q ~ m n /2, and that the phase transition line runs to that point directly, 
the curvature of the transition line for iso-vector /x should be much larger than that for usual 
/x, since m n /2 <C ton/3. However, as we discussed above, /x D in Eq.(27) does not contribute 
to the shift of ft near /x — and the difference from the usual /x is only in /x , i.e. /x G — for 
the iso-vector case. Therefore the difference in the curvature might be small and the naive 
picture seems to be wrong. In practice, our result at small /x using the method in Ref. [4] 
supports that. Moreover, Kogut and Sinclair [31] showed that ft from chiral condensate 
measurements is fairly insensitive to /x for small /x by direct simulations with iso-vector /x. 

B. Fluctuation of the reweighting factor 

Next, we estimate the fluctuation of the reweighting factor. As we have seen above, the 
fluctuation of the reweighting factor is separated into the complex phase factor of e AF and 
the other part, and these are almost independent. Moreover, this implies the absolute value 
of e AF is important for the determination of ft. The amplitude of the fermionic part |e AF | 
and the gauge part e AG are strongly correlated, and then the variation of the total fluctuation 
of these parts in the parameter space is not simple. Because the total fluctuation is related to 
the applicability range of the reweighiting method, here, we compute the standard deviation 
of |e F |e G to estimate the fluctuation, and also discuss the relation to the phase transition 
line. The complex phase fluctuation e td will be discussed in the next section separately. 

Up to the leading order of (3 — ft and /x 2 , the square of the standard deviation is obtained 

by 

m\e F \e G )f) « ((AP 2 ) V + 2(AP 2 AP)/x 2 (/3 - ft) + ((AP) 2 >(/3 - ft) 2 + • • • . (33) 

Then, the line of constant fluctuation is an ellipse in this approximation. We show the 
contour lines in Fig. 3. The susceptibilities and the correlation of R 2 and P, ((AP 2 ) 2 ), 
((AP) 2 ) and (AP 2 AP), are summarized in Table II. The values at the phase transition 
point, ft = 3.6497(16), are computed by the reweighting method for the f3 direction using 
the data at four f3 points. Numbers in this figure are the squares of the standard deviation 
divided by A^ site . We also denote the lower and upper bounds of <9 2 ft/<9/x 2 = —1.1(4) by 
bold lines, which are obtained by the measurement of the chiral susceptibility [4]. We find 



3 However, as /x increases, the difference between Xq an d 4xiv becomes sizeable [13], which might 
be related to only Xq being expected to have a singularity at the critical endpoint [29]. 
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that there exists a direction along which the increase of the fluctuation is relatively small, 
and this direction is roughly parallel to the phase transition line. Because we expect that 
physics is similar along the transition line, if we consider that le^e 6 * is the important part 
for the calculation of f3 c , this result is quite reasonable. 

As well as in the (m, (3) plane, the fluctuation of the reweighting factor is small along the 
phase transition line in the (fi, /3) plane, and the reweighting method seems to be efficient 
to trace out the phase transition line. This must be a reason why the phase transition line 
can be determined for rather large \x in Ref. [3]. However, in this discussion, we omitted 
the complex phase fluctuation, and the phase fluctuation is the most important factor for 
the sign problem. As we will discuss in the next section, the value of \i for which the sign 
problem arises depends strongly on the lattice size. The sign problem is not very severe for 
small lattices such as 4 4 , 6 3 x 4 and 8 3 x 4 lattices employed in Ref. [3], which is also an 
important reason for their successful calculation. 

Imaginary chemical potential In Fig. 3, we show also the region for fi 2 < 0, i.e. imag- 
inary fx. de Forcrand and Philipsen [10] computed d 2 f3 c /dji 2 performing simulations with 
imaginary /i, assuming that f3 c is an even function in fx and analyticity in that region. (Also 
in Ref. [11] for Nf — 4.) The f3 c {fx) for imaginary fx shifts the opposite direction to that for 
real fx as fx increases, but the absolute value of the second derivative, d 2 f3 c /dfx 2 1 is expected 
to be the same. Here, we confirm whether the results of \d 2 j3 c /d[x 2 \ obtained by real and 
imaginary fx are consistent or not by the method in Ref. [4]. We replace fx by ifx or —i/i and 
reanalyze for imaginary fx. In Ref. [4] , the reweighting factor has been obtained in the form 
of the Taylor expansion in fx up to 0(fx 2 ), and the replacement is easy. We determined f3 c by 
the peak position of the chiral susceptibility, using the data at itlq = 0.1 in Ref. [4]. The re- 
sults of \/3 c (fJ>) — f3 c (0)\ are shown in Fig. 4. Errors by the truncation of the Taylor expansion 
terms are 0(fx 4 ). The solid line is the result for real fi. The results of fi — > ifx and fx — > —ifx 
are dashed and dot-dashed lines respectively. The slope at fx — is ~-(d 2 /3 c /d[i 2 )/2. We 
find that these results of the slope for real and imaginary fx are consistent. It has been also 
discussed for measurements of spatial correlation lengths to confirm the reliability of the 
analytic continuation from imaginary chemical potential [32]. 

VII. COMPLEX PHASE FLUCTUATION 

Finally, it is worth to discuss the complex phase fluctuation in order to know the region of 
applicability of generic reweighting approaches. If the reweighting factor in Eq.(l) changes 
sign frequently due to the complex phase of the quark determinant, then both numerator 
and denominator of Eq.(l) become very small in comparison with the statistical error. Of 
course, the complex phase starts from zero at fx — but grows as fi increases. It is important 
to establish at which value of fi the sign problem becomes severe. 

As discussed in the previous section, the phase can be expressed using the odd terms of 
the Taylor expansion of In det M. The complex phase is 

9 = lm{R 1 fx + R 3 fx 3 + R 5 fx 5 + •••)• (34) 

The explicit expression for Rx and R 3 are given in Eqs. (24) and (26). 

Because the sign of the real part of the complex phase changes at 6 = ir/2, the sign 
problem occurs when the typical magnitude of 9 becomes larger than ir/2. We use the 
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point at which the magnitude of the phase reaches the value tt/2 as a simple criterion to 
estimate the parameter range in which reweighting methods will start to face serious sign 
problems. If the sign problem arises at small fj,, which is expected to happen for a large 
lattice, the first term in Eq.(34) is most important. Then we can estimate the applicability 
range by evaluating the fluctuation of R±. Moreover we expect naively that the magnitude 
of tr[M -1 (<9M/<9/i) • • •] is proportional to N site , therefore the value of /i at which the sign 
problem arises decreases roughly in inverse proportion to the number of site iV S it e . Also, 
the situation is different on lattices of moderate size. In Ref. [15], it is shown that the 
first term in Eq.(34) is dominant for /i — 0.1 and 0.2 but the higher order terms cannot be 
neglected for /x^;0.3, by calculating the complex phase without the approximation by the 
Taylor expansion. If the higher order terms are not negligible, the volume dependence is not 
simple. E.g. in the case that the term of 0(fi 3 ) plays an important role for the determination 
of the applicability range of the reweighting method, the applicability range is expected to 
decrease in proportion to N~^J 3 , and similarly N^J 5 for the case that the 0(/i 5 ) term is 
important. 

We consider the leading term and the next leading term of the complex phase. The 
expectation value of 9 must be zero at /j, — because the partition function is real. Although 
the average of the phase is zero, its fluctuations remain important. We investigate the 

standard deviation of 9 up to 0(/i 3 ), STD(#) = \J (9 2 ) — (6) 2 , using configurations generated 
on a 16 3 x 4 lattice for the study of Ref. [13], and also the standard deviations of Im(i?i) 
and Im(i? 3 ) are computed. 

The random noise method is used to calculate 9 for each configuration. Then, the value 
of 9 contains an error due to the finite number of noise vectors N Q . To reduce this error, we 
treat the calculation of (9 2 ) more carefully. Since the noise sets for the calculation of the 
two 9 in the product must be independent, we subtract the contributions from using the 
same noise vector for each factor. Details are given in the Appendix of Ref. [4] . By using 
this method, we can make the iV n -dependence of (9 2 ) much smaller than that by the naive 
calculation from rather small N n , hence it may be closer to the N n = oo limit. We took 
N n = 50 for this calculation. 

In Fig. 5, we plot the standard deviations of the 0(/i) term, Im(i?x), and the 0(/i 3 ) term, 
Im(i? 3 ), for Nf — 2, m — 0.1. The horizontal axis is temperature normalized by T c at /i = 
(To). The temperature scale is determined from the string tension data in Ref. [33] with the 
fit ansatz of Ref. [34] . The fluctuations of these terms are almost of the same magnitude and 
both of them are small in the high temperature phase, hence the sign problem is not serious 
in the high temperature phase. We, moreover, confirm that the 0(/i) term is dominant 
around /i = fi q a ~ 0.1, as suggested by Ref. [15], and the approximation by 0(/i) term 
for the discussion of the applicability range in Ref. [4] is valid for the 16 3 x 4 lattice. This 
suggests the applicability range decreases roughly in proportion to N^ e . However, in general, 
the magnitude of the fluctuation (standard deviation) of Ri/N site changes as a function of 
the volume, hence the detailed finite size analysis is necessary to investigate the volume 
dependence of the applicability range more precisely. 

Recently, analysis of the volume dependence of the applicability range has been reported 
in Ref. [35]. Their numerical result of the applicability range is in proportion to N si ^J 3 . This 
is much better than iV B 7 t g. Because their estimations are based on simulations on 6 3 x 4, 
8 3 x 4, 10 3 x 4 and 12 3 x 4 lattices, the applicability ranges are relatively large, and the 
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higher order terms in a should be important for large a. Therefore the result of the volume 
dependence, /i ~ N sit J , is reasonable for their lattice size, but may will change on large 
lattices. 

We also show the contour plot for STD(fl) = {tt/4, tt/2, 3vr/4, vr, 5vr/4, 3tt/2, 7tt/4, 2vr} in 
Fig. 6. The error of the contour is estimated by the jack knife method. At the interesting 
regime for the heavy-ion collisions: /i q /T c ps 0.1 for RHIC and /i q /T c m 0.5 for SPS, the 
fluctuation is smaller than tt/2 in the whole range of T. Therefore, the reweighting method 
seems to be applicable for the quantitative study for the heavy-ion collisions, which is an 
encouraging result. Also, we find that a point around T/T = 0.9 looks singular. Because 
we expect the fluctuation of the system diverges at a critical point, this might be related to 
the presence of a critical endpoint. The large fluctuation around T/Tq = 0.9, /j, q a > 0.5 is 
occurred by the 0(/i 3 ) term being large around there. This might be corresponding that, 
in general, the contribution from the higher order terms of the Taylor expansion becomes 
larger, as the critical point is approached, and the expansion series does not converge near 
the critical point. Plus sign and minus sign appear with almost equal probability, i.e. {e ld ) 
is almost zero, when the standard deviation of 9 is larger than n. The value of fi q /T = fiN T 
at which the standard deviation of 6 is tt around T c is of the order [i q jT ~ 0(1). However, 
we should notice that the complex phase, again, is very sensitive to the lattice size N site . 
For small lattice size, the sign problem is not severe and the reweighting method can be 
used for considerably large /x, however the applicability range of the reweighting method 
will be narrower for a lattice with a size larger than 16 3 x 4. Also the analysis of quark mass 
dependence must be important as a part of future investigations. 

VIII. CONCLUSIONS 

At present, the reweighting method is an important approach to the study of QCD at 
finite baryon density. We discussed the applicability range of the reweighting method with 
multi-parameters. The fluctuation of the reweighting factor during Monte-Carlo steps is a 
cause of the increase of the statistical error due to the reweighting, and the magnitude of 
the fluctuation determines the applicable range. 

For a simulation of SU(3) pure gauge theory at the first order phase transition point 
on an anisotropic lattice, the fluctuation is minimized along the phase transition line in the 
parameter space of an anisotropic lattice, if we assume the pressure in hot and cold phases to 
be balanced. This argument can also be applied in full QCD, if the theory has a first order 
phase transition, and the same relation between the fluctuation of the reweighting factor and 
the phase transition line is obtained. This suggests that the multi-parameter reweighting 
method is an efficient method for trace out the phase transition line in multi-parameter 
space. Moreover, the measurement of thermodynamic quantities on the phase transition 
line is important for the finite size scaling analysis to discuss the universality class. The 
reweighting method is useful for this purpose. 

We also measured the fluctuation of the reweighting factor in numerical simulations of 
2 flavor QCD for the cases of the reweighting in quark mass and chemical potential. There 
exists a direction of small fluctuation in the (m, j3) plane and it is roughly the same direction 
as that of the phase transition line. The fluctuation of the reweighting factor with respect to 
chemical potential can be separated into two parts: the complex phase factor of the fermion 
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part and, on the other hands, the absolute value of the fermion part and the gauge part. 
These two parts fluctuate almost independently during Monte-Carlo steps. This implies 
that the phase factor of the fermion part does not influence the shift of f3 c with increasing \i 
at small and also explains why the difference between the phase transition lines for usual 
chemical potential and iso- vector chemical potential is small at low density. If we neglect the 
complex phase factor, the increase of the fluctuation is also small along the phase transition 
line in the (fi 2 ,(3) plane, as well as in the (m, f3) plane. 

The value of \i for which the sign problem arises decreases as the lattice size A^ site in- 
creases, hence simulations at fi ^ are more difficult for larger lattices, even if the fluctuation 
of the absolute value of the reweighting factor is small along the phase transition line. The 
complex phase fluctuation is measured on a 16 3 x 4 lattice. The sign problem is not serious 
in the high temperature phase, but around the phase transition point it becomes serious 
gradually from this lattice size. For small lattices, the sign problem is not severe for the 
study at low density, and also for the 16 3 x 4 lattice, the applicability range of the reweight- 
ing method covers the interesting regime for heavy-ion collisions. Also, the behavior of the 
complex phase fluctuation around the transition point suggests a critical endpoint in the 
region of fi q /T c ~ 0(1). 

ACKNOWLEDGMENTS 

I wish to thank C.R. Allton, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, and 
C. Schmidt for fruitful discussion, helpful comments and allowing me to use the data in 
Refs. [4] and [13]. I also thank Y. Iwasaki, K. Kanaya and T. Yoshie for giving me the 
data in Ref. [19] and Ph. de Forcrand for useful comments on this manuscript. This work is 
supported by BMBF under grant No.06BI102 and PPARC grant PPA/G/S/1999/00026. 



15 



REFERENCES 



F. Karsch and E. Laermann, hep-lat/0305025 

S. Muroya, A. Nakamura, C. Nonaka and T. Takaishi, hep-lat/0306031. 

Z. Fodor and S. Katz, Phys. Lett. B534 (2002) 87; JHEP 0203 (2002) 014. 

C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, Ch. Schmidt 

and L. Scorzato, Phys. Rev. D66 (2002) 074507. 

Ch. Schmidt, C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch and E. Laer- 
mann, Nucl. Phys. B (Proc. Suppl.) 119 (2003) 517; hep-lat/0309116. 
Z. Fodor, S.D. Katz and K.K. Szabo, Phys. Lett. B568 (2003) 73. 
I.M. Barbour and A.J. Bell, Nucl. Phys. B372 (1992) 385; 

I.M. Barbour, S.E. Morrison, E.G. Klepfish, J.B. Kogut and M.-P. Lombardo, Phys. 
Rev. D56 (1997) 7063. 

PR. Crompton, Nucl. Phys. B619 (2001) 499; Nucl. Phys. B626 (2002) 228; hep- 
lat/0301001. 

A. Roberge and N. Weiss, Nucl. Phys. B275 (1986) 734. 

P. de Forcrand and O. Philipsen, Nucl. Phys. B642 (2002) 290; Nucl. Phys. B673 
(2003) 170. 

M. D'Elia and M.-P. Lombardo, Phys. Rev. D67 (2003) 014505. 

S. Choe et al. [QCDTARO collaboration], Phys. Rev. D65 (2002) 054501. 

C. R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and C. 
Schmidt, Phys. Rev. D68 (2003) 014507. 

R.V. Gavai and S. Gupta, Phys. Rev. D68 (2003) 034506. 

P. de Forcrand, S. Kim and T. Takaishi, Nucl. Phys. B (Proc. Suppl.) 119 (2003) 541. 
J. Ambjorn, K.N. Anagnostopoulos, J. Nishimura and J. J.M. Verbaarschot, JHEP 0210 
(2002) 062. 

A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 61 (1988) 2635; Phys. Rev. Lett. 
63 (1989) 1195. 

M. Fukugita, M. Okawa and A. Ukawa, Nucl. Phys. B337 (1990) 181. 
Y. Iwasaki et al, Phys. Rev. D46 (1992) 4657. 

S. Ejiri, Y. Iwasaki and K. Kanaya, Phys. Rev. D58, 094505; Nucl. Phys. B (Proc. 
Suppl.) 73 (1999) 411. 

F. Karsch, Nucl. Phys. B205 (1982) 285. 

G. Burgers, F. Karsch, A. Nakamura and I.O. Stamatescu, Nucl. Phys. B304 (1988) 
587. 

J. Engels, F. Karsch and T. Scheideler, Nucl. Phys. B564 (2000) 303. 
T.R. Klassen, Nucl. Phys. B533 (1998) 557. 

U.M. Heller, F. Karsch and B. Sturm, Phys. Rev. D60 (1999) 114502. 

S. Gottlieb, W. Liu, R.L. Renken, R.L. Sugar and D. Toussaint, Phys. Rev. D38 (1988) 

2888 

D. T.'Son and M.A. Stephanov, Phys. Rev. Lett. 86 (2001) 592. 
J.B. Kogut and D.K. Sinclair, Phys. Rev. D66 (2002) 034505. 
Y. Hatta and M.A. Stephanov, hep-ph/0302002. 

Y. Sasai, A. Nakamura and T. Takaishi, hep-lat/0310046. 

J.B. Kogut and D.K. Sinclair, Nucl. Phys. B (Proc. Suppl.) 119 (2003) 556; 
hep-lat/0309042. 

A. Hart, M. Laine and O. Philipsen, Phys. Lett. B505 (2001) 141. 

F. Karsch, E. Laermann and A. Peikert, Phys. Lett. B478 (2000) 447. 

C.R. Allton, hep-lat/9610016; Nucl. Phys. B (Proc. Suppl.) 53 (1997) 867. 

F. Csikor, G.I. Egri, Z. Fodor, S.D. Katz, K.K. Szabo, and A.I. Toth , hep-lat/0401016. 



16 



TABLES 



TABLE I. Correlation and susceptibilities of Qi and P. N sitc = 16 3 x 4. (3 C = 3.6492(22). 



p 


(AQ l AP)NZ 1 

\ ™ 1 / site 


((AQi) 2 )Nz\ 


((AP) 2 )Nzi 

W /I site 


3.640 


-1.44(6) 


1.29(6) 


2.67(7) 


3.645 


-1.80(13) 


1.69(13) 


2.99(14) 


3.650 


-1.70(6) 


1.58(6) 


2.89(7) 


3.655 


-1.76(14) 


1.65(13) 


2.92(15) 


3.660 


-1.60(6) 


1.49(6) 


2.77(7) 


3.665 


-1.36(13) 


1.19(12) 


2.58(14) 


3.670 


-1.52(7) 


1.41(8) 


2.68(8) 


3.6492 


-1.74(4) 


1.61(3) 


2.92(5) 



TABLE II. Correlations and susceptibilities among i?f,i?2 and P. N s \ te = 16 3 x 4. 
P c = 3.6497(16). 



p 


(A( J Rf/2)AP)iV s T t 1 c 


(A^AP)^ 


(A( J R?/2)AP 2 )iV & T t 1 c 




((AP) 2 )^ 


3.64 


0.006(29) 


0.312(33) 


0.034(10) 


0.216(16) 


2.62(10) 


3.65 


0.059(21) 


0.434(29) 


0.056(10) 


0.254(14) 


2.87(8) 


3.66 


0.055(15) 


0.410(26) 


0.022(5) 


0.231(11) 


2.75(8) 


3.67 


0.032(15) 


0.397(28) 


0.031(5) 


0.219(13) 


2.68(8) 


3.6497 


0.059(13) 


0.495(19) 


0.050(7) 


0.267(9) 


2.98(7) 
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FIGURES 




FIG. 1. Contour plot of the standard deviation of the reweighting factor and the phase transi- 
tion line in the ((3 a , (3 T ) plane. Bold line is the phase transition line, and the dashed lines denote its 
error. Values in this figure are the standard deviation divided by the mean value. The simulation 
point is ft = 5.6925. 




FIG. 2. Contour plot of the standard deviation of the reweighting factor in the (P, m) plane 
around P c , Pq = 3.6492, and mo = 0.1. Values in this figure are the square of the standard deviation 
divided by 7V s i te . Bold lines show the upper bound and lower bound of dp c /dm. 
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P-Po 

FIG. 3. Contour plot of the standard deviation of the reweighting factor in the (/3, /i 2 ) plane 
around (3 C , fio = 3.6497, and mo = 0.1. Bold lines show upper bound and lower bound of dfic/ddj, 2 ). 
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FIG. 4. |/3 c (/u) — P c (0)\ = |A/3 C | for real fi and imaginary fi. Solid line is the result for real \i. 
Dashed and dot-dashed lines are the results of \i — > i/j, and — > — ifi respectively. 
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